Polarons in the harmonic lattice. 
I. Standing polaron. 

V.N. Likhachev, T.Yu. Astakhova, G.A. Vinoghradov 
Emanuel Institute of Biochemical Physics, Russian Academy of Sciences, 
ul. Kosygina 4, Moscow 119991, Russia [*] 

Abstract 

We obtain analytical expressions for the large- and small-radius polarons on the one-dimensional 
lattice in the TBA approximation. The equations of motion for this model are treated classically for 
the oscillator subsystem, while a quantum description is used for the electron. The electron-phonon 
interaction is considered in the linear Su-Schrieffer-Heeger approximation. Good agreement be- 
tween analytical formulae and accurate numerical simulation is obtained. The dynamics of polaron 
formation from different initial conditions is considered. Some features of the wave function evo- 
lution, governed by the finite lattice length, are elucidated. 
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I. INTRODUCTION 



Charge and energy transfer is of utmost importance in both organic and inorganic nature. 
Studies of charge transport (CT) in organic systems were restricted till recently by research 
of amorphous state. CT in individual quasi-one- dimensional systems has attracted much 
attention. These oligo- and macromolecules have some advantages compared to traditional 
bulk systems because of less size, flexibility and power requirements. Biopolymers like 
natural and synthetic DNA, and proteins [l-4| are of particular interest. These systems can 
find an applications in solar cells, flexible TV displays, logical circuits etc. fl-Q. Molecular 



electronics can also be used in the nanotechnology lOHlq]. Photosynthesis, ATP hydrolysis, 



metabolism are examples where CT is vitally important 



Polarons are ubiquitous in materials where the electron-phonon coupling cannot be 
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nored. Polarons are thought to be responsible for the CT in biological systems 
including DNA. Different DNA models were considered to describe the polaron properties 
Dozens of papers present results on the effects of water and counterions surrounding 



DNA, on polaron properties in the presence of disorder, polaron hopping and drift in an 



applied electric field 
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Polarons in the harmonic lattice are thoroughly investigated in the resent paper. The 
electron-phonon interaction is described in terms of the linear Su-Schrieffer-Heeger approx- 
imation in the tight-binding modelling. Analytical expressions for small and large polarons 
are derived. Results coincide with numerical simulation with high accuracy. The dynam- 
ics of polaron formation from different initial conditions is considered. An influence of the 
lattice finite size on the wave function evolution is also analyzed. 



II. THE SMALL-RADIUS POLARON IN THE HARMONIC LATTICE 



A. Numerical simulation 



The one-dimensional harmonic lattice with one extra charge carrier is investigated. We 
consider an electron for the definiteness though the arguing is valid for the holes. The lattice 
consists of N particles. The "particle" represents one base of DNA in the coarse-grained 
model of DNA. The TBA approximation is used to describe the electron evolution on the 
lattice. 
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The hamiltonian of this system is the sum of the classical hamiltonian and the energy of 
electron-phonon interaction: 



H 



m 
~2~ 



N 



N-l 



(1) 



i=i 



i=i 



where Xi is the deviation of i-th particle from the equilibrium, \& = ipi, i/j 2 , ■ ■ ■ is the electron 
wave function, m and k are the particle mass and the lattice rigidity, respectively. In the 
TBA approximation the electron hamiltonian is the tridiagonal matrix: 

/ e x t x . . . \ 

tx e 2 t 2 ... 

t 2 e 3 t 3 ... 

H" = t 3 e 4 ... (2) 



... ejv-i t N -x 
\ . . . tjv_i e N } 

where t% is the hopping integral and e\ is the electron on-site energy. In the homogeneous 
lattice all diagonal elements are equal and can be set to zero, = 0. The hopping integral 
is written in the Su-Schrieffer-Heeger form {27!. I40I]: 



U = —[v Q - a(x i+ i - x^] 



(3) 



This model mimics the synthetic DNA or_peptide consisting of the regular sequence of 
identical bases or peptide groups (see e.g. [4l|). 

It is convenient to consider hamiltonian ([1]) in the dimensionless quantities. Let us choose 



the following units: vq is the energy unit, \Jmk 1 is the time unit and y/v^k 1 is the length 
unit. Then the dimensionless hamiltonian is 

^ N N-l N-l 

i=l i=l i=l 

where the same notation is preserved for the dimensionless variables; = (xj+i — Xj) is the 
relative displacement of the neighboring particles. The dimensionless parameter a/y/kvo — > 
a is the free parameter. Below we study the one- dimensional polaron solutions of the 
hamiltonian fl3j). 
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We find the solution for the standing stationary polaron when all velocities in (J3J) Xi = 0. 
The minimum of energy (j3J) with respect to all lattice variables {q{\ and the electronic wave 
function {ipi} determines the polaron. 

If Eq. (J3J) is minimized with respect to {ipi} then \l/ is the eigenfunction of the electronic 
hamiltonian H e with energy s. It means that two procedures, - energy minimization and 
solving the stationary Schrodinger equation, are equivalent. 

The electron wave function oscillates as ipie t£t / h , where {tpi} are real amplitudes of elec- 
tron wave function. Note, the time oscillations are eliminated in (j3J) because of complex 
conjugation. Then the expression for the energy (J3]) is 

^ N-l N-l 

H = -Y,ql-^Y. {1 - • ( 5 ) 

i=i i=i 

The minimization (jSJ) (i.e. dH/dqi = 0) gives 

% = -2aipiipi +1 . (6) 
Substitution of (jSJ) into d^J) results in: 

N-l N-l 

H = -2Y, AA+i ~ 2a 2 J](^+i) 2 - (7) 

i=J i=l 

An analogous approach was used in [2l| . 

Expression (J7|) was minimized numerically preserving the norm 

J2i^i = 1- If the lattice 

length iV is much larger then the polaron radius then the polaron solution does not depend 
on N. The numerical steepest descent method was used to find the polaron profile on the 
lattice with iV = 101. The result is shown in Fig. UK 

This polaron is centered on the particle. Such polarons will be called the site polarons 
hereafter. The polaron center can be located on the center of the bond between two nearest 
particles. Such polarons are called intersite polarons. The energy is somewhat higher in the 
latter case (see Fig. [lb). 

Small energy difference of the site and intersite polarons can be treated as the barrier 
separating two spatially different polaron states with equal energies. The typical value for 
the energy unit in DNA is Vo = 0.3 eV. Then the energy difference between site and intersite 
polarons is as small as < 0.0001 eV. This small energy difference points to the possibility of 
ballistically travelling polarons. This finding can also be the starting point for contraction 
of translationally-invariant states. 
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FIG. 1: Results of the numerical calculations of the polaron. a) Wave functions for three values of 
parameter a (shown in insert), b) Site (solid line) and intersite (dashed line) polarons for a 2 = 2. 
Energies of these polarons are shown in insert to the figure. Polarons are centered at the site i = 
(—50 < i < 50) in both figures. 



The total lattice contraction AL = g i; and using flH]) one finds 

N-l 

AL = -2a^2if>iip i+1 . (8) 
i=i 

For large-width polaron ipi « ipi+i an d taking into account the norm of the wave function 
(i.e. = 1) AL 2a . 



B. Analytical expressions for the small radius polaron 

For defmiteness we consider the site polarons hereafter. The wave function of polaron is 
symmetrical relative to its center. Consequently only one-half of a polaron (say, right half) 
can be analyzed. This problem is considered on the semi-infinite lattice and the new site 
numeration is used: i — 0,1,2, ... , oo. The amplitudes of the wave function are denoted as 
ipo for the central polaron site and ipi, i = 1,2,... for other sites. Then the expression for 
energy ((Tj) is transformed to 

oo oo 

H = -4 ^ Mi+i - 4« 2 J2(^i+if ( 9 ) 

i=0 i=0 
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FIG. 2: a) Exponential decay of the wave function for different values of a. N = 201. b) The 
dependence of parameter g vs. a. 

with the rescaled norm of the wave function 



Numerical minimization of energy ((9]) with respect to all ipi shows that the wave function 
decreases exponentially with the growth of i when i > 1 (see Fig. Efe>), he. ipi cx exp(— i/R) 
and R is the polaron radius. 

It is convenient to introduce parameter g related to the radius R by 



This relation follows from the exponential decay of the wave function amplitude: ipi oc g l at 
i > 1. In other words, g is the common ratio of the geometrical progression: g = ipi+i/ipi. 
Fig. [2b shows the dependence of g vs. parameter a of electron-phonon interaction. One can 
see that the value of g approaches unity for a — > and the polaron radius increases when 
parameter a decreases. 

The analytical solution can be constructed employing the asymptotic behavior of the 
wave function at large i. Our goal is to express the wave function and parameter a through 



It is convenient to introduce the Lagrange multiplier /i to find the energy minimum of 
(Q. It allows to take into account the condition for preserving the norm of wave function 
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(10) 



i=l 
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(|T0|) . Then the minimization of the functional 

/ oo 

# = tf + ^U 2 + 2^ 2 -l) (12) 



i=i 



with respect to ip , ipi, ip 2 , ■ ■ ■ gives the following system of equations 



(13) 



2^o 

+ ^i-i + 2a 2 ipi (tp 2 +1 + ijj 2 _ x ) = ntpu i > 1 
The Lagrange multiplier /x is expressed through g as 

A* ~ 9 + ~ • (14) 
9 

This relation follows from the second equation (TLB!) which becomes linear for i > 1: + 
= fiipi. The solution of this equation is ipi oc g l . It can be also verified that fi = — e 
where e is the electronic contribution to the polaron energy. 

The construction of the approximate analytical polaron solution will be performed in 
the following manner. We suppose that the exponential decay for the solution is a good 
approximation starting from some site number i . It means that the cubic terms in (fl3|) can 
be neglected. 

At the zeroth approximation the exponential decay starts from ipo (iq = 0) and ipi = ipo9 l - 
This is very rough approximation and we consider two next approximations. 



1. Polaron solution in the first approximation 

The exponential decay of the wave function starts from zq = 1 in the first approximation. 
And amplitudes are: i[) , ?pi,ip2 = V'i 9i ^3 = V'l • • •• Then only two equations for i[) and 
ipi remain in f fT3|) ; 

+ 2a 2 ifj ifjl = ^-ifj 

2 (15) 
^\9 + + 2a 2 ^i (V'o + 9 2 ipf) = A^i • 

An expression for the norm of the wave function is 

< + 2^, = 1 (16) 
and the second term in (Tl6l) is calculated as a sum of terms of infinite geometrical progression. 
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FIG. 3: a) Polaron energy vs. parameter a of electron-phonon interaction in the first and second 
approximations, b) The polaron radius R vs. parameter a. Results of numerical minimization are 
shown in empty circles. N = 101. 



As a result we get three algebraic equations (1151) and (TTb[) for three variables ipQ , ipi and 
a. The exact solution is too cumbersome and is not given. There exists a limiting case 
when term (gipi) 2 in (|T5l) is neglected. This approximation corresponds to the accuracy of 
the first approximation. The answer in this case is very simple: 
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(17) 



2 ' 2 

An expression for a is obtained by substitution ( TT71) and ( THl) into the first of equations ( |T5l) . 
The polaron energy, expressed through ipo,ipi and a, is 



E 



1-9 4 



(18) 



The comparison of the results obtained in the first approximation with the "exact" (nu- 
merical) results is shown in Fig. |3j There is some discrepancy and the next approximation 
is required. 



2. Polaron solution in the second approximation 

The exponential decay of amplitudes in the second approximation starts from i = 2. Am- 
plitudes of wave function are: ip , ipi, ip 2 , ^3 = ^2 9, ^4 = ^2 9 2 , ■ ■ ■ and form the geometrical 



8 



progression starting from the third term. Three equations 

' A* 

i/ji + 2a 2 ip ipl = 2 

from system (|T3|) is enough to solve the problem in the second approximation. The normal- 
izing condition for the wave function is 

^ 2 + 2^1 + = 1 . (20) 

Algebraic equations f|T9|) and fl20|) are solved numerically. The solution for the wave 
function is a series. First three terms (ipo, ipi, 1P2) are obtained from the solution of these 
equations. Next terms > 3) are members of the geometrical progression with the 

common ratio g. 

The comparison of "exact" (numerical minimization) and results obtained in the second 
approximation is shown in Figs. |3£i and [3b. One can see that the coincidence is very good. 



III. THE LARGE RADIUS POLARON 



The parameter a of electron-phonon interaction determines the radius of polaron: polaron 
has large radius if a is small. In this case the common ratio g in the geometrical progression 
is close to unity. We shall find the solution for the large-radius polaron in the continuous 
approximation on the infinite lattice (—00 < i < 00). 

The starting point for the evaluating the solution is the same as in the case of small-radius 
polaron. Let us write this equation once more for convenience: 

^+1 + Vi-i + 2a 2 ^ (</> 2 +1 + ipti) = M • (21) 

g is slightly less then unity for the large-radius polaron, and the Lagrange multiplier \x 
(fi = g + l/g) is somewhat larger then 2. Lets introduce small parameter 5 2 = /i — 2 and 
make some transformations in Eq. (l2Tj) . Firstly, subtract 2^ from both sides of Eq. |21j 

^ + A-i - 2il>t) + 2^-Vi W+i + Ci) = & • (22) 
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Next, (ipi + i + Vf-i) can be substituted by 2ipf because the difference between the amplitudes 
of the wave function on neighboring sites is very small for the large-radius polaron. As a 
result we get 

+ 4>i-i - 2^i) + = V>< ■ (23) 

And, finally, we introduce new variables: yi = — z—if>i- As a result the following equation is 

o 

obtained: 

^{y i+1 + y< _! - 2y<) + 2y? = y< . (24) 

First term in the left hand side of ( I24jl is the difference form of the second derivative. 
Applying the continuous approximation the following equation is obtained: 

y" = y + 2y 3 . (25) 

Its solution is 

1 6 1 

V = — rrr ° r *l>(x) = — j= — rr— . (26) 
cosh (a;) ay 2 cosh J 

The normalization of the wave function gives 5 = a 2 , then the final expression for the 
large-radius polaron in discrete form is 

a 1 , , 

^ = 7Icosh(a 2 *)' { } 



Very similar solution is obtained for a somewhat analogous problem in [23]. 

The comparison of "exact" results (numerical minimization) and the analytical expression 
(|27|) is shown in Fig. H]for two values of parameter a. The agreement is very good and it is 
the better the larger is polaron radius. In one dimension case there is a transition from the 
small to the large polaron and a ~ 0.6 separates these two solutions. 



IV. THE DYNAMICS OF POLARON FORMATION 

Two cases of large and small polarons were considered above. These solutions correspond 
to the energy minimum in the stationary state. And now we consider the evolution of the 
dynamical and electronic degrees of freedom giving the polaron as the final state. 

The following numerical experiment was performed. The wave function with some profile 
was chosen as the initial electronic condition. Initial conditions for the dynamical variables: 
lattice at rest, i.e. Xi = 0, ii = 0. 
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FIG. 4: The dependence of the amplitudes of wave function vs. the site number i according to 
(|27p . Dashed line - for a = 0.4, solid line - for a = 0.2. Circles - accurate results of numerical 
minimization. 



The system evolution is determined by the following system 



where qi 



oti =(xi-i + x i+1 - 2xi) + a[(ipiip* +1 + ip*ipi+i) - (ipi-iip* + ipi-^i)] 

• i (28) 
ipi =~ [(1 - ai&M-i + (1 - aqi)i/)i + i] 
h 

. — - Yl 1 777, . — 

(xi+i — Xi) and h is the dimensionless Planck's constant: — * > h. If DNA 

v o V « 



parameters are chosen 20| then 1/h ~ 80. 

If the wave function is initially located on one site (^(t = 0) oc ^ ), or has the constant 
value on the lattice (il>i(t = 0) oc A^ -1 / 2 ), then the wave function spreads over lattice very 
rapidly and polarons are not formed. The reason is that the characteristic electronic time is 
very small (~ h ~ 10 -2 ) compared to the dynamical time scale (~ 1). And the dynamical 
degrees of freedom have no enough time to follow the wave function. 

If the initial wave function is chosen as a eigenfunction of unperturbed lattice, then 
the lattice "captures" the oscillating electron and the joint electron-lattice evolution re- 
sults in polaron formation. The trial complex wave function is chosen in the form 
${t = 0) = | v/— -sin 



N + 1 



Til 



N+1 



. It corresponds to the ground-state wave func- 



tion of the unperturbed lattice. The evolution of the wave function is shown in Fig. |5j The 
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FIG. 5: Snapshots of the wave function evolution at different time instances shown in insert. 
Initial wave function ipi(t = 0) oc (sm[iri/(N + 1)],0). a 2 = 2, N = 101. At t = 14 the polaron is 
completely formed. 



polaron is formed during few oscillation periods of the lattice. 

Eqs ( l28|) correspond to the microcanonical ensemble where total energy is preserved. The 
energy gain due to the polaron formation is compensated by the vibrational excitation of 
a lattice. Because of the dynamical interaction with the lattice, polaron fluctuates, slightly 
changes its form and position, but stays very stable. If frictional forces at the lattices ends 
—7x1 and — 7Xat are introduced then all vibrations decay and only polaron survives on the 
lattice. 

If the initial wave function is the eigenfunction of the first excited state ip}(t = 0) = 

I then two polarons are formed but with norm of wave function 



N + l 



sin 



2tu 



N+l 



= 0.5 for each of them. 

The problem of charge transfer from a donor to the DNA double strand is very complex 
quantum-dynamical problem. An accurate solution of of time dependent Schrodinger equa- 
tion for coupled state of excited donor and a lattice is necessary. We did not touch this 
problem here and limited ourself by the model representation of the initial wave functions. 
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V. WAVE FUNCTION EVOLUTION ON THE FINITE-SIZED LATTICE 



Polarons on the infinite lattice were considered in the previous sections. But if the 
particular applications of the suggested approach are of interest (e.g. charge transfer in 
DNA) then the lattice should have the finite length, i.e. N < 100. And there appearers few 
peculiarities which differ the polaron dynamics on the infinite lattice from the lattice of final 
length. Few examples are considered below. 

First of all, the intriguing property of the evolution equations (|28|) should be mentioned. 
At some initial conditions the interaction of an electron with the lattice is the "one-way": 
There exists an influence of the lattice on the electron through electron-phonon interaction, 
but the back "polaronic" action of an electron on the lattice is absent. In other words, 
additional forces acting on particles from the electron subsystem are zeroes. In the simplest 
case it is realized if the wave function is initially localized on one (arbitrary) site. Let the 
initial wave function for definiteness is real and totally located on the first site: ip\{t = 0) = 1. 
From the second set of equations (1281) it follows that at the certain time step the wave 
function on odd sites is real, and on even sites is imaginary. At the next time step the 
situation inverts: the wave function on odd sites is imaginary, and on even sites is real. 
Then the products V'iV'i+i * s always zero ( see (|28|) ) . Hence, the lattice always stays to be 
harmonic. This conclusion is also valid if the initial wave function is organized in such a 
manner that it is real on all even sites and imaginary on all odd sites (or vise verse). 

Second interesting case is the wave function evolution on the inhomogeneous lattice. Let 
one (arbitrary) lattice site i$ has the TBA parameters which differs from others. These 
parameters are ^ for the electron on-site interaction, and a ^ a for the electron- 
phonon interaction fl3]). Initial wave function can be chosen arbitrary. If the wave function 
for the homogeneous lattice (e = 0,a = a) never forms the localized state - polaron, but 
for the lattice with defect the case is opposite: the polaron is formed after some induction 
period. The result is shown in Fig. Though the norm of the polaron wave function is 
less then unity, polaron stays a stable quasiparticle for a very long time. 

The polaron formation can be observed as follows: The maximal value of the wave func- 
tion amplitude is monitored and this maximum is averaged over some time interval. The 
formation of the stable maximum with max|?/>| = 0.5 — 0.6 is the direct evidence of the 
polaron formation (see Fig. |6b). Polaron is formed practically instantaneously at 5 = 1, 
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FIG. 6: a) Snapshot of the polaron wave function at time t = 134. a 2 = 2, N = 30,e25 = 
1-5, a25 = 0.1. b) Maximal value of the wave function (in absolute value) vs. time. Solid line 
for a = 0.1, dots for 5 = 1. Every point is the result of averagin on the time interval t = 1.25. 
N = 30, a 2 = 2. 



and after induction period t « 100 for a = 0.1. 

There exists one more interesting case of the wave function evolution at some special 
choice of parameters. If parameter e = (absence of the diagonal disorder) and only 
one hopping integral differs from the others, then it is the case when the lattice does not 
influence the electron evolution (as discussed above). The usual "hopping" dynamics on 
the homogeneous lattice will be observed. Let only first site differs from the others by the 
hopping constant t\ ^ t (see (j2J)) and the wave function is localized on the site i — 1 at 
the initial time t = 0. Naturally, no polaron is formed as the lattice does not influence the 
electron dynamics. But there exists the phenomenon of the periodic returning of the wave 
function to the first site. This phenomenon is related with the final lattice size and is shown 
in Fig. [7J Time is measured in "electronic" time units ht which are related to dynamical 
time units by the dimensionless Planck's constant h = 0.0125 (see above). 

This problem (in very simple formulation) reminds experiments on the photoinduced 
charge transfer from intercalated metallocomplex to the DNA chain 42J], and the value t\ 
reflects the interaction of the ligand with the DNA bases. 

Some interesting effects are experimentally observed: "chemistry at a distance" 4^ and 
"piug-poug e.ect r o„ transfer Q Our resu.ts can help to elucidate some experiments 
results which has no satisfactory interpretation yet [2jj. Particularly the detailed analysis of 
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FIG. 7: Absolute value of the amplitude of the wave function on the first site. N = 30, ti = 0.1. 

the "electronic ping-pong" will be given in our next paper where moving polarons will be 
considered. 

The polaron behavior and wave function evolution were considered at zero temperature till 
now. The kinetic energy of thermal fluctuation can be rather high at ambient temperatures, 
and the polaronic time of life can be limited due to interactions with these fluctuations. 
Below we consider the evolution of the wave function at finite temperatures. 

The following problem was solved numerically to find out the role of temperature. The 
stationary polaron was formed at the lattice center, as described in the previous section. 
Then the lattice was thermalized from its ends (Langevine random forces with friction) 
to achieve the necessary temperature. If the temperature is high enough, the polaron is 
destroyed rather soon. We monitor the maximal value of the squared absolute value of 
the wave function amplitude on the lattice, computing its averaged value over time interval 
t = 50. In the beginning this value is close to ~ 0.4, what approximately corresponds to the 
squared maximum of amplitude for the standing polaron. Then this maximum falls down 
to ~ 0.1 at some time r. This time corresponds to the polaron life time. When temperature 
decreases, r increases. Life time r vs. temperature is shown in Fig. [8] 

Of course, there no any firm basis to make a conclusion that the found temperature 
To is the true critical temperature. Sooner this temperature should be determined by the 
lattice length N and the parameter a of electron-phonon interaction. It seems reasonable 
that the "critical" temperature should decrease with the growth of the lattice length N and 
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FIG. 8: Polaron time of life r vs. the temperature. Temperature is measured in units of energy 
(vq = 0.3 eV for DNA). Solid line is the hyperbolic fitting lg(r) = a + b/(T — To), where a, b are 
some parameters, and To is the "critical" temperature. Dashed line shows the asymptotic value 
To ~ 0.016. Error bars correspond averaging over 20 random trajectories. 

decreasing parameter a. 

VI. CONCLUSIONS 

Analytical expressions for the large- and small-radius polarons are derived in the present 
paper. The polaron radius is determined by the parameter a of the electron-phonon inter- 
action in the SSH model: polaron is small at a > 0.6. The solution for the wave function of 
small-radius polaron represents a series where the first three terms are the solution of the 
system of algebraic equations and the next terms are geometrical progression. This form of 
solution is explained by the exponential decay of the wave function with increasing distance 
from the polaron center: ipi+i/ipi = g = const. For the large-radius polaron an analytical 
expression is derived in the continuous approximation. This solution reminds the solutions 
of some nonlinear equations (nonlinear Schrodinger, modified Korteweg - de Vries). Good 
agreement between analytical solutions and numerical simulation is obtained. 

The considered simple model can help to explain existing experiments on the long range 
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coherent charge transfer in DNA. 
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